function [dependent_variable,regressor_mat,time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs,reg_matrices]...
            =create_regression_matrices_no_NaN_indicator(dependent_variable,regressor_mat,data_array_for_regression_stacked_by_variable,...
                    pos,country_indicator_names_mapping,time_fixed_effects_dummy,country_fixed_effects_dummy,indicator_mat)
% function [dependent_variable,regressor_mat,time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs,reg_matrices]...
%            =create_regression_matrices_no_NaN_indicator(dependent_variable,regressor_mat,data_array_for_regression_stacked_by_variable,...
%                    pos,country_indicator_names_mapping,time_fixed_effects_dummy,country_fixed_effects_dummy,indicator_mat)
% Create regressor matrix by dropping all time-points that have no non-NaN observations
% We need at least one country with full set of regressors
% 
% Inputs:
%   - dependent_variable                                [T by n_countries] dependent variable matrix
%   - regressor_mat                                     [T by n_countries by nvars] regressor matrix before removing NaNs and adding fixed effects
%   - data_array_for_regression_stacked_by_variable     [T by n_countries by nvars] data matrix including fixed effects
%   - pos                                               [structure] containing variable position
%   - time_fixed_effects_dummy                          [boolean]   indicator whether time fixed effects are requested
%   - country_fixed_effects_dummy                       [boolean]   indicator whether country fixed effects are requested
%   - indicator_mat                                     [T by n_countries] indicator matrix
% Outputs:
%   - dependent_variable                                [T_non_NaN by n_countries] dependent variable matrix after removing NaNs
%   - regressor_mat                                     [T_non_NaN by n_countries by nvars] regressor matrix after removing NaNs and adding fixed effects
%   - time_indices_without_only_NaN_obs                 [T_non_NaN by 1] indices of non-NaN variables
%

% Copyright (C) 2019-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% For a copy of the GNU General Public License,
% see <http://www.gnu.org/licenses/>.

if nargin<8
    error('Not enough input arguments')
else
    if ~isequal(size(dependent_variable),size(indicator_mat))
        error('Matrix dimensions do not fit.')
    end
end
drop_warning_dummy = 0;
time_dim=1;
country_dim=2;
var_dim=3;

%% Stack matrices
n_dependent_variables=size(dependent_variable,3);
stacked_matrix=cat(3,regressor_mat,dependent_variable,indicator_mat);
%get NaN-entries
NaN_entries=isnan(stacked_matrix);
%find country-time observations with only incomplete set of variables
country_time_NaN_obs=any(NaN_entries,var_dim);
%find timepoints without country with  full set of variables
timepoints_without_full_set=all(country_time_NaN_obs,country_dim);
time_indices_without_only_NaN_obs=find(~all(country_time_NaN_obs,country_dim));

%select unproblematic time-points
stacked_matrix_without_NaN_timepoints=stacked_matrix(~timepoints_without_full_set,:,:);

%% get countries to be dropped; NB: cannot affect time points, because only countries without full set of regressors are dropped
NaN_entries=isnan(stacked_matrix_without_NaN_timepoints);
country_time_NaN_obs=any(NaN_entries,var_dim);
country_without_full_set=all(country_time_NaN_obs,time_dim);
country_indices_without_only_NaN_obs=find(~country_without_full_set');

if any(country_without_full_set) && drop_warning_dummy
    fprintf('\n')
    country_index=find(country_without_full_set);    
    fprintf('I am dropping %s because there are only NaN observations in the sample\n',country_indicator_names_mapping{data_array_for_regression_stacked_by_variable(1,find(country_without_full_set),pos.country_ident)})
    fprintf('\n')
end

stacked_matrix_final=stacked_matrix_without_NaN_timepoints(:,country_indices_without_only_NaN_obs,:);


regressor_mat=stacked_matrix_final(:,:,1:end-1-n_dependent_variables);
dependent_variable=stacked_matrix_final(:,:,end-1-n_dependent_variables+1:end-1);

n_countries=size(dependent_variable,2);

%% add fixed effects
data_array_for_regression_country_Fixed_Effects=[];
if country_fixed_effects_dummy
    for ii=1:length(country_indices_without_only_NaN_obs)
        temp=data_array_for_regression_stacked_by_variable(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs,pos.(['Country_FE_' num2str(country_indices_without_only_NaN_obs(ii))]))...
            .*indicator_mat(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs);
        % cannot identify time-fixed effect based one country observation.
        if sum(sum(temp,'omitnan'),'omitnan')>1
            data_array_for_regression_country_Fixed_Effects=cat(3,data_array_for_regression_country_Fixed_Effects,temp);
        end
        temp=data_array_for_regression_stacked_by_variable(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs,pos.(['Country_FE_' num2str(country_indices_without_only_NaN_obs(ii))]))...
            .*(1-indicator_mat(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs));
        if sum(sum(temp,'omitnan'),'omitnan')>1
            data_array_for_regression_country_Fixed_Effects=cat(3,data_array_for_regression_country_Fixed_Effects,temp);
        end
    end
end
data_array_for_regression_time_Fixed_Effects=[];

if time_fixed_effects_dummy    
    if country_fixed_effects_dummy %drop one time FE to avoid dummy variable trap
        end_range=length(time_indices_without_only_NaN_obs)-1;
    else
        end_range=length(time_indices_without_only_NaN_obs);
    end
    for ii=1:end_range
        temp=data_array_for_regression_stacked_by_variable(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs,pos.(['Time_FE_' num2str(time_indices_without_only_NaN_obs(ii))]))...
            .*indicator_mat(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs);
        NaN_index_country_time=any(isnan(stacked_matrix_final),3);
        temp(NaN_index_country_time)=NaN;
        % cannot identify time-fixed effect based one country observation.
        if sum(sum(temp,'omitnan'),'omitnan')>1
            data_array_for_regression_time_Fixed_Effects=cat(3,data_array_for_regression_time_Fixed_Effects,temp);
        end
        temp=data_array_for_regression_stacked_by_variable(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs,pos.(['Time_FE_' num2str(time_indices_without_only_NaN_obs(ii))]))...
            .*(1-indicator_mat(time_indices_without_only_NaN_obs,country_indices_without_only_NaN_obs));
        NaN_index_country_time=any(isnan(stacked_matrix_final),3);
        temp(NaN_index_country_time)=NaN;
        if sum(sum(temp,'omitnan'),'omitnan')>1
            data_array_for_regression_time_Fixed_Effects=cat(3,data_array_for_regression_time_Fixed_Effects,temp);
        end
    end
end
if ~country_fixed_effects_dummy && ~time_fixed_effects_dummy %add constant to regression
    data_array_for_regression_time_Fixed_Effects=ones(size(regressor_mat,1),n_countries);
end

reg_matrices.regressor_mat=regressor_mat;
reg_matrices.data_array_for_regression_country_Fixed_Effects=data_array_for_regression_country_Fixed_Effects;
reg_matrices.data_array_for_regression_time_Fixed_Effects=data_array_for_regression_time_Fixed_Effects;
reg_matrices.dependent_variable=dependent_variable;

regressor_mat=cat(3,regressor_mat,data_array_for_regression_country_Fixed_Effects,data_array_for_regression_time_Fixed_Effects);
